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Quantitative Non-Destructive Testing May 1 , 1984 - April 30, 1985 

Executive Summary 

The work undertaken during this period included two primary efforts. 
The first is a continuation of theoretical development from the previous 
year of models and data analyses for NDE using the OPTITHIRMS system, 
which involves heat injection with a laser and observation of the resulting 
thermal pattern with an infrared imaging system. The second is an 
investigation into the use of the thermoelastic effect as an effectove tool 
for NDE. As in the past, the effort is aimed towards NDE techniques 
applicable to composite materials in structural applications. 

The theoretical development started from the ideal cases described^ 
earlier and produced several models of temperature patterns over several 
geometries and material types. Agreement between model data and 
temperature observations was obtained. A model study with one of these 
models investigated some fundamental difficulties with the proposed 
method (the primitive equation method) for obtaining diffusivity values in 
plates of finite thickness and supplied guidelines for avoiding these 
difficulties. A wide range of computing speeds was found among the 
various models, with a one-dimensional model based on Laplace’s integral 
solution being both very fast and very accurate. 

Thermoelastic experiments were successfully performed on a stainless 
steel plate with a hole in the center, and some relevant system 
characteristics were obtained. A literature review found related activity 
in several other places including foreign countries. One of these activities 
has led to a commercially available device for making thermoelastic 
stress measurements, and a demonstration was held. A possible 
application to NDE in previously strained materials was suggested by the 
data. The situation for graphite-epoxy composites was found to be 
complex both experimentally and theoretically compared to that for 
metals, and directions for further study are suggested. 
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1 . 0 Introduction 

This report discusses the work done under NASA Cooperative Agreement 
NCCI-50 with the Physics Department of the College of William and Mary for the 
contract year from May 1984 through April 1985. As none of the work has been 
contained in progress reports for portions of this period, this report will 
cover the topics included in a complete manner. 

The scope of the reported work includes two primary efforts. The first 
is theoretical support of the OPTIcal THermal ^Infra-Red Measurement System.^. 
(OPTITHIRMS) , which is being developed at NASA Langley for non-destructive 
examination of structural materials using infra-red radiative techniques and 
extensive computer-based data analysis of the resulting measurements. The 
emphasis of this work is on graphite-epoxy composite materials. The second 
effort involves the use of the thermoelastic effect to examine stresses in 
materials undergoing cyclic, controlled loading, with the goal of identifying 
flawed and weakened portions of structures prior to failure. This effort is 
also focused on applicability to graphite-epoxy structures. It has both 
theoretical and experimental elements. 

2.0 Thermal Models and Diffusivity Estimates 

The present configuration of the OPTITHIRMS produces a heated area on the 
surface of a sample by irradiating it with a carefully scanned beam from an 
infrared laser. The heat pattern thus produced is radiatively detected with a 
scanning infrared imager, and the image data is digitized and stored for later 
analysis. To date, this system has been applied to estimates of the 
diffusivity of a material from the time evolution of the resulting temperature 
pattern. 

The models developed to date have been in support of this effort . In the 
previous year, the models consisted of analytical models for simple 
geometries, and the emphasis was on qualitative issues, such as sensitivity as 
a function of dimensionality and appropriate non-dimensionalizations of the 
conduction equation. It was found that if time is scaled to the pulse time of 
heat input, a natural spatial scale is generated by combining the time scale 
with the material diffusivity. In infinite or semi-infinite media, the 
heating curves for point sources then reduce to products of physical 
constants, such as emissivity, heat capacity and pulse power, and universal 
functions of scaled time and space. These functions were derived and 
illustrated for configurations with dimensionalities of 1,2 and 3. The effort 
during the current year was to apply these idealized heating curves to 
experimentally observed data . 

The computer techniques used to formulate these models were more 
sophisticated than those for the previous year. In the previous year, the 
formulas were obtained in closed form in terms of standard transcendental 
mathematical functions. Numerical computations were used to evaluate 
polynomial approximations to these functions at desired points in time and 
space. In the present year, the source functions became distributed, and 
additional numerical techniques were used to sum partial temperatures over the 
source at each desired response point for a given value of time. The models 



thus constructed are termed numerical integration models. 

In the area of data analysis, the first estimate of diffusivity of a 
material was obtained by image analysis based of evaluation of terms directly 
from the differential equation for heat conduction. The method of calculating 
diffusivity directly from estimates of terms in the governing equation is 
called the primitive equation method. As experimental use of this method was 
found to be difficult, a model study was undertaken to examine the application 
of the primitive equation method to plates of finite thickness, and guidelines 
for controllable experimental parameters required to obtain accurate 
diffusivity values were obtained. 

2.1 Line scan models 

The motivation for the line scan model was to synthesize a temperature 
field from theoretical considerations which would be directly comparable to 
radiometric data obtained from the OPTITHIRMS system. For this purpose, 
motivated by considerations of theoretical simplicity, experimental 
simplicity, and closeness of correspondence, the situation chosen was the 
surface temperature distribution on the face of a thick sample of stainless 
steel on which a short, straight line of heat was injected with the laser. 

The surface of the stainless steel was coated with a high-emissivity coating 
(emissivity > .99) to boost signals, but the coating was thin enough that 
viewed temperature corresponded closely with surface temperature of the 
sample. Experimental durations were chosen to be short enough that the sample 
appeared semi-infinite, with reflections from the sides and back being well 
below the limits of detection. The short line was chosen to boost the maximum 
temperature attained in the heated area, again to boost the signal level. The 
model assumed a Gaussian distribution for the beam profile. Parameters for 
the model included beamwidth, line length, line position and orientation, 
pulse duration, input power, sample emissivity, image scale, and image time. 
Material properties required were thermal diffusivity, conductivity, and 
volumetric heat capacity. Calculations were obtainable over both the heating 
and cooling phases of the experiment. 

The model, based on the three-dimensional point source equation, divided 
the image surface into its individual pixels (squares of size A) and treated 
each of these as a point source located at the pixel center. The continuous 
local variables, x and y, thus became referenced in integral blocks of size A 
with integer variables I and J respectively. The temperature at the (I, J) 'th 
block thus was represented by the sum of partial temperatures from each (K,L) 
block with non-negligible heating intensity, represented by the intensity at 
the center of the block J(x K ,y L)j with x K = A(K - 1/2) and y L = A(L - 1/2) . 

With these definitions, the distance between a source point and an observation 
point is given as r IJ>KL = A[(I-K) 2 - (J-L) 2 ] 1/2 . The temperature, is thus 
formally represented as: 
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in which the source is pulsed between times (t) zero and t off and U is the 
Heaviside step function. In local coordinates, with the zero in the center of 
the source line, which extends in the y direction for a total distance of L, 
and has a scale width of w. the source intensity function is 
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This formulation is appropriate except when I=K and J=L, the response pixel 
and the source pixel being the same. In this case, the value of r IJKL becomes 

zero, and the corresponding term becomes infinite. The artifice of assigning 
some nominal value to this case and letting the pixel size become very small 
is not satisfactory because convergence is very slow and, particularly for 
short times, local temperature is primarily due to local heating. Thus, the 
model must be altered in this case, and the most reasonable alteration is to 
abandon the notion of a point source and spread the integration somehow over 
the local pixel. In doing this, because only the local pixel is considered 
and because the pixel size is generally much smaller than the scale length for 
the Gaussian source (A/w « 1), it is also appropriate to use the asymptotic 
time- independent formula for temperature distribution within the pixel and to 
consider the intensity within the pixel to have a constant value equal to the 
value at the pixel center, I. Here, use is made of the concept, valid only in 
three-dimensional geometries, that the temperature pattern asymptotically 
approaches a time-independent shape rather than increasing without limit under 
constant heat input. With a constant source intensity limited to a finite 
region, the formulation of Lax [1977] shows the temperature at an observation 
point to be proportional to the average of the reciprocal distance between 
observation and source averaged over the source region. This average may be 
denoted <l/r>. The temperature contribution calculated in this manner is not 
suitable as a representative contribution to the total heating pattern from 
the local pixel because it depends on the particular position of the 
observation point within the pixel. To obtain a representative contribution, 
it is necesary to average this temperature over all of the observation points 
within the pixel, an operation denoted by «l/r». For a square pixel with 
sides of length A, the formula for «l/r» is 
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The family of integrals in this expression, which can be called 1^, can be 
expressed in terms of a single member 1^ with the scaling expression 1^ *= A 3 ! 1 . 
This integral was evaluated in an approximate manner by reducing it to a sum 
of integrals over subvolumes of side 1/N. Only a small number of these 
integrals contain the singularity, and these can be expressed in terms of 1^ ' 
If this is done and the remainder of the integrals are approximately evaluated 
as the product of their 4-dimensional volume and the central value of the 
integrand, the formula 
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results, the primes on the sums denoting that the singular terms (those with 
i=k and j=l) have been removed. After further manipulation to consolidate the 
number of summations from 4 to 2, the expression was evaluated for an 
increasing series of values for N up to N=500, and an approximate value of I* 

= 2.971 was obtained as the asymptotic value approached for large N. This 
value was then used to provide the representative value for the singular pixel 
during the period of heating. A different approximation was necessary within, 
this pixel to account for partial temperature contributions during cooling (t 
> t off) , for the use of the asymptotic value in the difference formula which 
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applies following the pulse results in zero contribution to the local pixel 
from heat applied at the same pixel. To make this approximation, the error 
function terms in the point source expression were expanded in Taylor Series 
to second order in their local times. As expected, the first order terms 
summed to zero, leaving a contribution at second order. 

This resulted in the formula: 
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for the partial temperature at a pixel identical with the source pixel at a 
time subsequent to the end of the laser pulse. Inspection of this formula 
reveals that at the moment following the end of the pulse, the partial 
temperature from this formulation again becomes infinite. Thus, there is a 
"blind time" in the model for a short period following the pulse. This blind 
time corresponds to the time prior to the establishment of the asymptotic 
time-independent temperature within the pixel. The blind time was estimated 
as the time for the contribution at the source pixel following the pulse to 
equal the contribution immediately before the end of the pulse. This equality 
leads to an estimate of the blind time of the same order as A 2 /100k. For 
stainless steel and the .0115 cm pixel size considered in the verification 
effort, the blind time was on the order of .0001 s. With a frame time of .03 
s, the blind time for the verification was considered unimportant for all 
frames following the first after the visible laser pulse. 

Model formula development thus resulted in three separated cases, which 
provided a complete description of the temperature evolution in a descrete 
grid, with the exception of a short time, the "blind time", following the end 
of the pulse of heat input. The formulas describe, respectively, the 
contribution for all pixels for which source and response locations are not 
identical, the case of identical source and response pixels during the time of 
the pulse, and the case of identical source and response pixels during the 
time following the pulse. An additional formula describes the duration of the 
blind time following the pulse . 

Model construction, using these formulas, was first done on a grid of 480 
pixels vertically and 512 pixels horizontally, to correspond to the digitized 
image data. In order to construct an evaluation algorithm, one further issue 
needed resolution. The integrals in the Green's function evaluation have 
infinite limits, a clearly impossible situation for numerical evaluation of 
sums. It becomes necessary, then, to describe a finite set of pixels for 
summation by identifying a locale within which all of the source pixels which 
have a significant effect on a given response pixel are located. At the same 
time, a criterion for "having a significant effect" needs to be defined, for 
if the effect is too small, the underflow errors will occur during the 
calculation, and computer time will be devoted to calculating insignificant 
contributions. For the heat equation, this task is simplified somewhat, as 
the temperature resulting from a given localized heat input decreases 
monotonically with distance in all directions from the source. In the 
homogeneous material covered by this model, the requirement may be met by 
defining a "radius of significance" around each source point beyond which the 
heat from the source has no significant effect. Sources are only considered 
if they are closer to the response point under consideration than the radius 
of significance. This criterion would be sufficient if all source points had 
equal intensity, but this is clearly not the case for the laser line scan 
considered by the model. There is, in addition, a region around the source 
outside of which no significant heat is added to the sample. From the 
iso-intensity curves surrounding the source function, this region is seen to 
have the shape of a capsule with sides parallel to the laser scan and 
semi-circular "end caps". For a given response pixel, contributions to the 



temperature rise are computed only for source pixels located within both the 
capsule around the laser line and circle with the radius of significance 
around the response point. When the two regions have no intersection, the 
calculation is terminated. In addition, the calculation is limited to those 
response pixels lying within a "window of interest", chosen as the area within 
the digitized image. For the model, the radius of significance was defined in 
terms of a user-chosen permitted fractional error. The region chosen was the. 
smallest set of pixels which completely covered the circle outside of which 
the total fractional response was equal to the permitted error fraction. The 
source capsule was chosen to be the smallest set of pixels covering the 
iso-intensity contour outside of which the total fractional incident power is 
equal to the permitted fractional error. 

The model which evaluated these equations subject to these constraints 
for a given user-chosen laser line and material was written in FORTRAN 77 as a 
program called LINETEMP with its associated subroutines. The program is 
resident on the MCIS VAX and on backup tapes. A further discription is given 
in Appendix 1 . The output of the model is a single file containing the 
calculated temperature field as an array of real numbers, with associated 
date, time and array size values. This file was written to be commpatible 
with a conversion program to put it into SADIE (The image analysis and display 
package resident of the MCIS VAX) format for display and further analysis. 

A significant property of the program was its running time on the MCIS 
VAX 750 in timesharing mode. The verification run involved a line pattern with 
a length of about 240 pixels, half the height of the image window, and a 
width of about 100 pixels, about 20% of the image window width. The evaluation 
of this image required about 50 hours of computer time, clearly too much for 
practical application. On analysis of the operation time, over 95% was 
associated with forming the sums, with the other 5% or so used in evaluating 
the expressions in the formulas. The time-consuming element in summation was 
"page faults", the disk accesses required in a "virtual memory" computer for 
obtaining values separated too far in address space. This separation distance 
is directly a function of the size of the arrays used in the model. To 
ameliorate this situation, the model was reduced in resolution by a factor of 
4, to produce a region of interest 128 pixels by 120 pixels. This reduction, 
besides providing substantial increases in speed, is compatible with the 
manufacturer's specifications for the number of independent data points 
sampled over an image by the infrared imager. The expected running time was 
estimated to be proportional to the number of term in the sums. For a given 
image, the number of terms is proportional to the fourth power of the 
resolution. Reducing the resolution by a factor of 4 was thus expected to 
reduce the run time by a factor of 256 to about 12 minutes. The actual 
runtime attained was about 6 minutes, the further reduction being attributed 
to a smaller number of disk accesses required with the smaller (l/16th the 
original size) arrays. Further reduction could reasonably be expected by 
reconfiguring the runtime environment of the model to provide more main memory 
during operation. 

For model verification, data were taken from a short line scan of a block 
of stainless steel, and the model was run for nominal laser and stainless 
steel parameters corresponding to the experimental conditionss. The data were 
limited in value by the absorbtion of heat in the stainless steel sample, 
resulting in a relatively low signal to noise ratio. In order to obtain a 
reduction in noise, a band was averaged across the transverse part of the 
image to include 60 original points in each of the averages. The temperature 
scale was not available in the data, so the range of averaged temperatures was 
set equal to the temperature range in the model. The central positions were 
likewise estimated from the image data. The resulting model verification data 
are shown in Figure 2.1. In this figure, the temperature profile of the data 
is seen to compare closely with that of the model. 
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Figure 2.1. Comparison of data and line scan model for stainless steel 
thick sample. Comparison is between a transverse profile of the image data 
comparably extracted from both model and data images. Variations between 
the model and the data are within the limits expected from noise in the 
observations . 

The thin plate situation has several advantages over the thick sample 
case. The temperature distribution from a given source continues to increase 
with irradiation time, rather than reaching an asymptotic limit, so the 
temperature contrast in an image can be increased simply by increasing pulse 
length. Also, the distribution of heat within the sample equilibrates 
throughout the thickness of the sample in a time of order L 2 /k following the 
cessation of thermal input, where L is the sample thickness and k is its 
diffusivity. Thus, the temperature distribution throughout the sample 
subsequent to this time is completely determined by the surface temperature 
distribution. This permits specification of the thermal state of an object 
from sets of measurements by the same thermal imager over time, eliminating 
the requirement for consistent measurements of laser pulse characteristics and 
material emissivity and the implied associated requirement for 
intercalibration of diverse sensors. Further simplification is attained by 
placing a laser stripe on the thin plate, for then the situation is reduced to 
one spatial dimension if the sample is homogeneous.. Finally, the reduction in 
dimensionality permits models to describe the thermal field to within a given 
resolution and acccuracy with greatly decreased runtimes. As will be shown is 
subsequent reports, these advantages have permitted considerable progress 
using linear heat patterns on thin plates. 

During the period covered by this report, a model called TLINE was 
written to use an initial one-dimensional set of temperature data as its input 
and produce sets of temperature predictions at times subsequent to the initial 
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data time as its output in the case of a cooling thin plate following a 
line-scanned input. The basis for the model is Laplace’s solution of the 
one-dimensional heat equation for an infinite solid, adapted from Carslaw and 
Jaeger [1959, eq 2.2(1)] as: 

oo 
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where x is position in the one-dimensional space, k is diffusivity of the 
material, and t is the time since the initial conditions were observed. 

In evaluating this integral, the spatial coordinate is divided into intervals 
of common length A, so that the distances (x-x') all become integral multiples 
of A, jA. Then all of the arguments of the exponentials become j 2 powers of a 
common factor, F x = exp (-A 2 / (4kt) ) and the sum strongly resembles a 

convolution sum. Finally, an iteration scheme was used by which all of the 
factors can be constructed from products of previous factors, the result being 
that the exponential function need only be called once for each calculation 
time. The single set of exponential factors and the iterative calculation of 
these factors with only a single evaluation of the exponential function both 
save considerable computer time. Further time could be saved if the sums were 
cut off after the terms become smaller than some cut-off value, but this 
strategem has not yet been implemented in this model. With real data, values 
of the initial conditions do not extend to infinity, of course. Still, some 
value needs to be supplied to the formula beyond the ends of the line at long 
times to approximate these unknown values, or else the ends of the line will 
start to exhibit artifacts. In TLINE, the option exists for specifying a 
constant value to be assumed for the initial temperature beyond the ends of 
the data. If a value of zero is supplied, the program uses the end data 
points for the constant value. 

In tests against a Gaussian curve, TLINE reproduced the analytic values 
to better than 1 part in 10 6 under the main peak. In the edges, artifacts of 
up to .005 degrees have been found out of values on the order of 10 degrees. 
This numerical accuracy is presently greater than the experimental data by 
several orders of magnitude, so further refinement seems inappropriate at this 
time. Running time for a single time prediction of 256 points on the VAX 750 
in MCIS is 1.7 seconds. Because it requires only one parameter, diffusivity, 
and makes few assumptions about the input data, program TLINE should prove 
useful whenever temperature data are to be predicted in one-dimensional 
situations from previously existing temperature data. Tests have shown that 
the algorithm is not suited for back-prediction, as small perturbations grow 
very quickly to dominate the solution. 

2.2 Point source temperature distribution 

Work from the previous year showed that the time evolution of a point 
source in a three-dimensional, semi-infinite medium approaches an asymptotic 
spatial distribution varying as the inverse of the distance from the source. 
With a real beam, having its intensity distributed over some finite area, this 
implies that the central temperature in the heated zone will approach some 
finite value, as shown by Lax [1977] . For a given beam shape, the attained 
temperature is thus dependent only on the total intensity, which can be 
controlled very precisely with laser sources and presicion attenuators. This 
precise control of the maximum attainable temperature has important 
implications in the materials processing industry. It can also be used in the 
context of NDE for purposes of calibration as well as for thick sample 



scanning. To exploit this feature fully, it is desired to obtain an 
expression for the radial dependence of the asymptotic temperature profile 
rather than simply its peak value. To this end, the starting point is the 
point source temperature distribution in a semi-infinite three-dimensional 
medium, given by Carslaw and Jaeger [1959, eq 10.4(2)]: 

T= ^Kr erto (?4kr) 

As time approaches infinity, the complementary error function approaches 
unity, reducing the functional dependence of the temperature from a point 
source to only the input power and the distance from source to observation 
point. The calculation of the asymptotic temperature distribution from a 
distributed source is accomplished by weighting the inverse distance from each 
source point to a given response point by the power density at the source 
point and integrating over the entire source pattern. For a circular Gaussian 
source pattern with a halfwidth of w. Lax [1977] expresses the result as 

T - T max N(R,Z,oo) 

where R and Z are the radial distance and depth normalized by the source 
halfwidth: R = r/w, Z = z/w. The third parameter in the N function is related 
to the absorbtion depth of the incoming radiation. Its infinite value 
corresponds to surface absorbtion, the usual assumption in NDE of structural 
materials. With this notation, T max = P/(2VrcKw) incorporates the direct 
effects of the power and the beam width, K being thermal conductivity of the 
material, while N is a function of the non-dimensional position. The image of 
the heated spot is due only to surface temperature, so our interest is in 

N (R, 0,°°) . From Lax [1977], this function is the ratio of two integrals, 

fj 0 R)exp(- X 2 /4) d X 
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The numerator is evaluated using the tables of Gradshteyn and Rhzhik [1980, eq 
6.618], while the denominator is a common form. On evaluation of N and 
substitution of the expression for T max the asymptotic variation of the 

temperature from a Gaussian spot is given as 

T= 2^r ex P ( -2^ )I '> ( 2^> 

where Iq is the zero order modified Bessel function of the first kind. This 

particularly simple form for the limiting temperature distribution from a 
Gaussian heat source in a thick material has not previously been reported. 

Its numerical evaluation is relatively fast, for the modified Bessel function 
is well described by a 6 term polynomial in the square of the argument for 
small arguments, while for large arguments, the product of the modified Bessel 
Function and the exponential function are well described by a 9 term 
polynomial [Abramowitz and Stegun, 1965, eq 9.8.1, 9.8.2].0 

2.3 Point source model and Diffusivity Estimates 



The development of a "point source model" was undertaken in order to 
examine the application of the primitive equation technique with the thin 
plate assumption for estimating diffusivity in plates of finite thickness. 

The approach taken was a Green's Function approach similar to that for the 
line scan. Considerable computation was saved, however, by calculating the 
resulting temperature only along a single radius, the material under 
consideration being homogeneous and the source radially symmetrical. The 
sample thickness is taken into consideration by constructing a series of image 
sources, each of which adds to the temperature at a point within its region of 
influence. Because of this method of construction, the point source model is 
applicable both to plates of finite thickness and to very thick plates. It is 
coded as program POINTEMP, and has as output a time series of temperatures 
between two chosen times and at a chosen time interval between two chosen 
radial distances at a chosen radial interval. 

The method of diffusivity estimation used during the first part of the 
contract year was direct evaluation of ratios of temperature derivatives in 
the heat equation itself. This method is known as the primitive equation 
method, to identify it among several methods which have since been introduced. 
From the experimental side of the effort, the question of how long the source 
should be left on to obtain good diffusivity estimates was posed. The answer 
to this question divides into two parts. The first part concerns the 
temperature field itself and addresses estimation issues associasted with 
uncertainty in signal levels and souces of noise. The second part of the 
answer concerns applicability of the formulation itself: In a plate of 

finite thickness, under what conditions does the ratio of observed temperature 
derivatives approximate the material diffusivity? It is the second part of 
the answer which motivated the model development. 

The first approach to this answer was to construct an analytical 
expression for the temperature on a finite-thickness plate from a series of 
point sources, the real source and the set of image sources associated with 
the three-dimensional distribution of temperature in an infinite medium which, 
taken as a sum, satisfy the adiabatic boundary conditions used for the flat 
plate. The only difference between the total solution to the heat equation 
and that obtained by evaluating terms visible on the surface only is that the 
depth dependence of the temperature is neglected in the evaluation of the 
surface evaluation. The ratio of surface-evaluated to true diffusivity is the 
ratio of the entire Laplacian to the 2-dimensional (surface) Laplacian. The 
actual calculation for the error was that of the ratio of the difference of 
these two to the entire Laplacian. The expectation with this approach was 
that a set of lines could be constructed in the position-time domain each of 
which corresponded to a given error, say 5%,1%,.1% ..., and that at long times 
following the laser pulse, as the temperature field equilibrated through the 
plate thickness,, the error would asymptotically approach zero. The result was 
different from this expectation. It showed that for any chosen distance from 
the source, the error never fell below some finite value, and that this value 
increased monotonically on approaching the source point. Stated differently, 
it showed that there was no finite region where diffusivity could be 
determined with any desired accuracy simply by waiting longer for the 
temperature to equilibrate through the plate thickness. 

The way to reconcile the non-convergence of the point source analysis 
with the seemingly trivial identity in the two-dimensional case seemed to be 
to distribute the energy in the point source over some finite area on the 
plate and rerun the image solution with sources distributed over the image 
planes rather than point sources. Because the difficulty was considered to be 
associated with the singular nature of the point source, distributed sources 
with continuous derivatives were expected to cure the problem. A source 
distribution with Gaussian form would not only satisfy this criterion, but 
would also correspond to the best available estimate of the distribution in 
the actual laser beam used in the experiments. Thus the error estimates 



obtained from the numerical results would be directly applicable to the 
experimental situation. A search of the literature revealed no analytic way 
of adding solutions to the heat equation at layers of images with Gaussian 
distribution, so a numerical model was called for. The resulting model was 
developed in a program named POINTEMP . 

During the time that the theoretical work had been going through these 
phases without producing a satisfactory answer to the question of how long to 
make the laser pulse, the experimental part of the work moved to the use of 
line scans rather than point scans. An exhaustive survey of parameter space 
with POINTEMP was considered relatively unimportant compared with work on 
models with line sources, so only a qualitative survey was done with POINTEMP 
and its diffusivity estimator, CALCK. The diffusivity estimator uses a 
numerical approximation to the time derivative and the Laplacian from the 
radial data produced by POINTEMP. The only notable feature of CALCK is that 
the time derivative estimator is fashioned to use the same data values with 
the same weights as the Laplacian estimator. Thus no errors in calculated 
dif fusivities are attributable to an offset in the locus of the two estimates. 
The estimates employ six points in radius-time space and may be considered as 
centered at the central of three spatial points and halfway between the two 
time points. Figure 2.2 shows the results of a single run of POINTEMP and 
CALCK for stainless steel. 
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Figure 2.2. Model results showing regions of high and low accuracy for 
difffusivity calculations carried out with the primitive equation analysis on 
the surface of a homogeneous plate of finite thickness. The plate, with a 
diffusivity of .04 cm 2 /s and a thickness of 2 mm represents the stainless 
steel laboratory samples. For this run, beam radius was 1 cm, pulse duration 
was 2 seconds, and power was 8 Watts. Peak temperature under these conditions 
at the end of the laser pulse reached 10.3 degrees C. 


This run, although for only a single case, illustrates the behavior of the 
regions of both high and low accuracy. It is a plot of the accuracy of the 
diffusivity calculations as a function of both radial distance from the center 


of the source and of time, and it includes all of the time for which good 
accuracy is obtained from the primitive equation method for diffusivity 
calculations. During the time that the laser is on, the diffusivity 
calculation gives accurate results in the far field, in corroboration with the 
analytic calculation. Immediately following the laser pulse, there is no 
region of accurate calculation. This corresponds to the L 2 /k period during 
which the thermal pattern equilibrates through the thickness of the sample. . 
Following this time, there appear two regions in which the calculations 
achieve accuracy of better than 1% separated by a band of inaccuracy. These 
will be called the near field and far field regions. The far field region 
recedes away from the source as time passes, while the near field is present 
for only a limited time. In the case presented, the width of the laser beam, 

1 cm, is much greater tban the plate thickness, 2 mm. Other tests show that 
the width of the laser beam spot must be greater than the sample thickness for 
the near field to achieve high accuracy. 

The qualitative behavior of the calculation can now be seen for the 
relationship between the distributed source case, which shows a limited region 
of high accuracy, and the point source solution, which shows only a far field 
region of high accuracy. During the heating time, the horizontal distribution 
of heat is largely determined by the source characteristics . In the far 
field, the temperature pattern is close to that expected from a point source " 
of the same power, so the far field region of accuracy is established for the 
distributed case in the same way as for the point source. This region of high 
accuracy is likely to be of little use in NDE, because of small signal levels. 
Immediately following the source pulse, the heat flow is dominated by 
equilibration through the thickness of the plate, so the primitive equation 
method measures essentially irrelevant spatial derivatives, the result being 
low accuracy. In the time following this equilibration, which seems to be 
about half of the scale time (L 2 /k) , a high accuracy region is found to be 
located directly beneath the heat source. The radial limit of this region may 
be explained because, in an expanding heat spot, the center is cooling while 
the outer part of the pattern is heating. Thus, there is a location where the 
time derivative of temperature passes through zero. If the primitive equation 
calculation were perfect, the measured horizontal Laplacian would be zero at 
this same place. The numerator and denominator of the primitive equation thus 
both become zero in this region, so non-zero values from the calculation 
reflect entirely the small numerical inaccuracies in the model. With real 
data, they would reflect noise. At any rate, the calculation fails in a dead 
band around this zero time derivative region. As time increases, this dead 
band increases to the point where it eventually consumes the near field region 
of high accuracy. Subsequent to this, the accuracy of the solution resembles 
that found for the point source calculation. This corresponds to the concept 
that, at long times following a confined source pulse, the temperature from 
the pulse increasingly resembles that from a point impulse with equivalent 
total energy. In the case examined, the region of high accuracy persisted 
during a period between 0.5 and 1.5 times the scale time and within a radius 
equal to the acale redius of the source beam. 

Based on the foregoing analysis, the answer to the original question can 
be expressed as a rule-of-thumb. To obtain good diffusivity estimates with 
the primitive equation method, the laser beam should be broadened to cover an 
area with diameter substantially larger than the thickness of the sample. 

This broadning need not preserve the Gaussian figure of a focused laser beam, 
but should have a central maximum and no "doughnuts" in the pattern. It 
should be directed at the sample for long enough to produce, following 
equilibration through the thickness, a good signal-to-noise ratio. If a scale 
time is defined as the ratio of the square of the nominal sample thickness to 
its nominal diffusivity, good accuracy in diffusivity measurements can be 
expected in the period from one-half to three-halves of the scale time in the 
portion of the sample under the heating pattern. If a nominal value of 



diffusivity is not available, scale time can be estimated as the time for the 
original rapid temperature decrease following heating to occur under the 
central peak of the temperature pattern. 

From a mathematical perspective, there is a behavioral similarity between 
the time variation of the diffusivity calculation and the accuracy of the 
asymptotic expansion for the exponential integral with an increasing number of 
terms [Morse, P. M. and H. Feshbach, Sec. 4.6]. If one considers the terms 
which are significant in the POINTEMP program, each one corresponds to a 
small source of heat located within the heated region on one of a set of image 
planes. Taken together, the sources are located within a bounded region at 
any given time, the bounded region having a shape like a stack of pancakes. 

At short times, few of the images are included in the sums, corresponding to 
just a few pancakes, the stack having mostly the shape of a single pancake. 

As time increases, more and more pancakes are added to the stack, which 
eventually resembles a long cylinder. At very long times, the cylinder 
resembles a string, the original shape of the pancakes becoming a relatively 
unimportant aspect of the total shape. The string corresponds to the shape of 
the analagous source region for the point source approximation. The accurate 
solutions are found within this progression during the time when the bounded 
region appears to contain several pancakes, or image planes, but before the 
time when it resembles a stack of coins. Thus, the region of accuracy for the 
diffusivity estimate does approach the correct diffusivity in an asymptotic 
manner, the approach reversing following some finite number of terms. In this 
sense, the two-dimensional model for thermal flow on a thin plate is seen to 
be an asymptotic approximation to the more realistic finite plate calculation. 

2 . 4 Model Summary 

The emphasis during the contract year shifted from examination of 
idealized models directed towards elucidating fundamental properties of 
thermal flow to models which can be compared directly with OPTITHIRMS data. 

The models constructed were all of the numerical integration type, forming 
sums of responses from patterns of idealized sources. Modifications to the 
formulation of the numerical equations were made to accommodate the 
singularity at zero radius. Model running time and computer environment were 
found to be significant new issues. To illustrate, two models which were 
constructed to produce temperature distributions from line scanned laser 
sources were found to have running times of 51 hours and 1.7 seconds, 
respectively. Besides providing correspondence with the available 
experimental data, the models were used to examine the correspondence between 
the two-dimensional idealization of a thin plate and its more realistic 
three-dimensional counterpart with respect to the primitive equation method 
for determining diffusivity. 

3. Thermoelasticity 

The thermoelastic effect is the thermodynamically associated change in 
temperature of an elastic material subjected to stress under adiabatic 
conditions. It is frequently ignored in elastic analyses because it 
represents a small part of the total stored energy in most common fabrication 
materials, and the temperature variations are quite small, on the order of 
hundredths of a degree . In contrast to the heat produced by material 
yielding, thermoelasticity is a reversible effect and is generally expressed 
as cooling under tension, rather than heating. Because the level of detection 
with the OPTITHIRMS system can be greatly increased under cyclic loading and 
because the thermoelastic effect produces, in principle, images which carry 
information about the stress distribution in an object under elastic cyclic 
loading, an investigation was initiated under this grant into the possible 
utility of thermoelastic image data for nondestructive evaluation of 



graphite-epoxy composite structural members. In contrast with existing 
methods of examining stress distributions, thermoelastic variations detected 
radiometrically have the potential of providing image data, as opposed to 
collections of point data, with minimal disturbance of the object under test, 
for they sample only the thermal radiation which is emitted anyway by the 
object 

The work under this phase of the contract was performed in collaboration 
with Sidney Allison of the' Materials Characterization and Instrumentation 
Section of NASA Langley Research Center. He provided access to and operation 
of the MTS Stress-Strain machine in the section as well as a wealth of 
knowledge concerning mechanical testing techniques, sample clamping methods, 
and associated matters. 

3 . 1 Literature Review 

The first theory of thermoelasticity was formulated in 1855 by William 
Thomson (Lord Kelvin) , and experimental verification was obtained to better 
than 0.1 percent in 1915 by Compton and Webster [Oliver, et al, 1982]. In 
oscillating motions, Zener [1938] applied the thermal conduction arising from 
thermoelastic temperature variations to the problem of damping of vibrations . 
Biot [1956] combined the thermomechanical and conductive equations to show 
that, when material deformation associated with temperature changes in 
conducting materials is taken into account, it is entropy rather than 
temperature which satisfies the diffusion equation. The actual use of 
infrared radiometry to examine structural stress was first reported by Belgen 
[1967a, b] . Another report of experiments relating stress to thermoelastic 
cooling was given by Heyman and Chern [1982]. Nayroles, et al. [1981] noted 
the existence of thermoelastic oscillations in infrared data taken to observe 
the propagation of a fatigue crack, and remove the thermoelastic signal as 
part of a search for the region of plastic deformation producing heating at 
the crack tip. Bouc, et al. [1983] reported the use of digitization of the 
infrared images obtained from a thermal imager to produce an image of the 
trace of the stress tensor around a crack tip. Oliver, et al [1982] reported 
the use of an infrared thermoelastic device which builds up images from a 
series of point measurements obtained through synchronous demodulation of 
signals with a stress-correlated driving signal associated with a sample under 
cyclic stress in its elastic regime. They made several comparisons between 
experimental results and theoretical solutions, both analytic and finite 
element, for the cases of a hole in a flat plate and a cylinder under 
compression. This work has led to the introduction of a product, the SPATE 
(Stress Pattern Analysis by Thermal Emission) 8000, a device which makes the 
synchronous detection approach to thermoelastic stress analysis commercially 
available. 

The work reported above has been applied entirely to homogeneous 
materials, primarily metals. For this case and at temperatures well away from 
absolute zero, such as room temperatures commonly used in NDE, the temperature 
fluctuation under adiabatic conditions is directly proportional to the product 
of the ambient temperature and the sum of the principal stresses. The 
constant of proportionality is the ratio of the linear thermal expansion 
coefficient to the volumetric heat capacity at constant pressure, and the sign 
is chosen so that tension produces cooling. These relations are expressed as 

6 = -K m TO (3.1) 

with 

K M - C^/Cp. (3.2) 

Here 0 is the amplitude of the temperature fluctuation, o is the amplitude of 
the sum of the principal stresses, and T is the absolute value of the ambient 


temperature. Kfj i.3 called the thermoelastic coefficient, and it is a property 
characteristic of any given material, is the linear coefficient of thermal 

expansion, and Cp is the volumetric heat capacity at constant pressure. 

Belgen [1967b] notes that the actual temperature fluctuation is different at 
high stresses, which are frequently used to obtain large signal values in 
samples, because the pressure within the sample varies with the strain during 
each cycle, and so Cp is not the appropriate value, but rather an intermediate 
value between Cp and c v . 

It is not immediately clear how this formula should be altered for the 
case of non-isotropic materials, for the temperature fluctuations evidently 
come about through an interaction of the strains with the thermal expansion 
terms. In an anisotropic material, both are represented by second order 
tensors, and the appropriate interactions are not intuitively clear. 

3.2 Experimental Work 

The experimental work has been carried out in facilities made available 
through the Materials Characterization Information Section (MCIS) of NASA 
Langley Research Center. The particular experimental configuration consisted 
of an MTS fatigue stress machine, rated at 110 thousand pounds, the VAX 750 
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Figure 3.1. Experimental configuration for thermoelastic measurements. 


computer associated with MCIS, and an Inframetrics Model 525 Infrared Imager. 
The VAX is attached to a Grinnell image processing system with a video 
digitizer and to an LPA-11 laboratory interface, which is used as a 
synchronization pulse generator. The experimental configuration is shown in 
Figure 3.1. 

In operation, the imager generates the video synchronization signal, 
which cycles 30 times per second and i3 not synchronized with the remainder of 
the system. The MTS machine is set up with a single cycle ramped pulse 
initiated from an external synchronization pulse. The ramp times, hold time, 
and maximum load level are set with panel controls on the MTS controller. The 
VAX, under the program MTS_Control (Appendix 1) , initiates a data cycle by 
instructing the Grinnell digitizer to obtain a positive image and add it to 
its image buffer. After a user-chosen delay, the VAX instructs the Laboratory 





Interface to send the synchronization pulse to the MTS machine. The computer 
also sends, following the synchronization pulse by a user-chosen interval, 
another instruction to the image processor to digitize a negative image and 
add it to the image buffer The cycle is repeated after a settling time 
subsequent to the second image. The cycle timing diagram is shown in Figure 
3.2. In addition to running the data cycle, MTS_Control initializes the 



Figure 3.2. Thermoelastic cycle timing and control. Max load, ramp time 
(a) and hold time (b) are selected on the MTS Fatigue Loading Machine. 
The cycle starts with the instruction to acquire a positive image, which 
is actually obtained at the next imager frame cycle (2/60 sec) 
subsequent to the command. The MTS cycle is initiated after a 
synchronization time (c) following the cycle start. The negative image 
is obtained at the first frame time following a time interval (d) from 
the sync pulse, and the cycle is repeated after a settling time (e) 
following the negative image. Cycle time is the sum of (c) , (d) and (e) . 


Grinnell by setting its image buffer to an intermediate value and halts 
following a user-chosen number of cycles, saving the final image in a disk 
file with a user-chosen name. Later versions of the program also transfer the 
contents of the Grinnell image buffer to a summing buffer in the VAX following 
a small number of cycles, also user-chosen. The reason for the transfer is 
that the image buffer sums and stores quickly, in a frame time, but has 
limited dynamic range. The computer memory has very large dynamic range, but 
transfer from the Grinnell requires 9 seconds per image. The sub-sampling 
scheme is chosen to provide the large dynamic range of the computer buffer to 
the data while incurring a minimum time penalty from the data transfer 
operation. To increase the effective dynamic range at the expense of 
resolution, the program also shifts the incoming data down one bit following 
digitization and prior to summation. The low order bit from the 8 bit 
digitizer (1 part in 256) contains little information when the camera noise 
level is .15 degrees in 10 (3.8 parts in 256), so its loss is of little 
significance. 

On the advice of the technical staff of the Materials Division of NASA 
Langley Research Center, who are thoroughly versed instress distribution, the 
equipment was used to examine a sample consisting of a metal plate with a 
circular hole in the center, as a close approximation to a situation with 
plane stress distribution around a hole in an infinite plate subject to 
uniform unidirectional stress at large distances away from the hole. This 
case, besides applying to holes used for fastening sheet metal, is 
experimentally convenient because the sum of the principal stresses is 
described by a simple function which is localized around the hole and attains 



both positive and negative values 

<T - <t,£l " 2 £|Jcos (20) J (3.3) 

Here, R is the radius of the hole, the formula applying only outside this 
radius, 0 is the angle of a polar coordinate system with zero angle in the 
direction of far field stress, denoted as r is the radial coordinate for 

the polar system, and C is the sum of the principal stresses . 

The experiment to examine the system response compared with the ideal 
response was performed on a plate of Type 347 stainless steel 2.95 mm in 
thickness. Stainless steel was chosen because its relatively small value of 
thermal diffusivity reduced the evolution of the temperature field due to 
thermal conduction during the period of production and acquisition of 
experimental data. Thus, the adiabatic approximation was experimentally 
valid. The plate was 102 mm wide and 629 mm long with a hole 19.0 mm in 
diameter placed in its center. Its central portion was coated with a high 
emissivity coating to produce high-efficiency emission in the infra-red. The 
plate, which was self-aligning in the MTS machine owing to its width in the 
hydraulic grips, was cycled for 20 times at a maximum load of 5,000 pounds. 

The load cycle had ramp times of 50 ms and a hold time of 100 ms. The cycle 
repeat time was 990 ms. Image acquisition instructions were issued at cycle 
times zero and 200 ms. The images obtained were saved, and the experiment 
repeated 5 times. As four of the five realizations looked qualitatively 
similar in the resulting stress images, their sum, containing a total of 80 
stress image pairs, was used to form a stress image with random noise reduced 
by a factor of 8.9 over a single image pair. The fifth realization was 
rejected because of a qualitative judgement that image "banding" due to 
instrumental noise was too great. 

The image data were further processed to obtain the temperature deviation 
profiles along the principal axes of the plate, and the data are shown 
compared to theoretical values in Figure 3.3. The two pairs of data/model 
curves are taken from the image data along the two principal axes of the 
experiment, through the hole along the center line of the plate and across the 
plate perpendicular to the center line and through the center of the hole. To 
reduce noise, the data were first filtered by averaging over a 4 by 4 pixel 
area from the video image. The resulting average values each cover an area 
equal to about two independent points of infrared data, the translation scheme 
used to obtain video formatting being highly redundant. In forming the data 
lines, the samples are further averaged over four parallel lines nearest to 
the principal axes, the lateral variation being small in these regions. 

The theoretical values used in Figure 3.3 were taken from equation 3.3 
and modified using a filter to describe the sampling area of a single image. 
The filter was a running average (a "boxcar" filter) encompassing the points 
within a given segment length of the model. The segment length is chosen to 
match the slope obtained in passing the filter over the edge of the hole with 
that observed in the data, under the assumption that the radiation field at 
the edge of the hole approaches a step function as closely as possible. 

The fit of the experiment to the model in figure 3.3 is qualitatively 
good, indicating that the temperature variations observed in the data are due, 
as anticipated, to the thermoelastic effect. Of particular interest is that 
the area of net compression predicted from the theory corresponds to the 
observed temperature rise in the same region. The data thus exhibit 
temperature fluctuations in the same general pattern as those expected from 
thermoelastic theory and have signs indicating both heating and cooling, 
consistent with thermoelasticity and inconsistent with dissipative 
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Figure 3.3. Comparison between thermoelastic data and model. Fitting 
parameters include 67 vertical units for the temperature depression 
within the hole, the zero value; 85 vertical units for the temperature 
depression corresponding to unit stress; and .12 inches for the 
detector "footprint" size. Comparisons for both along-sample and 
across-sample axes are shown, with the observations showing the 
effects of instrumental noise . 

heat-producing effects, the other candidates for producing temperature 
contrast in the image. 

There are aspects of the data which invite further interpretation and 
analysis, however, as they are to a greater or lesser extent questionable or 
paradoxical. The stress pattern extending across the plate dips down at the 
edges of the plate (at plus and minus 2 inches) on both edges. This suggests 
some kind of edge effect which may be of interest or may reflect some 
experimental factor not presently included. The size of the boxcar window 
required for a close fit to the data came to about .12 inches, or 1/32 of the 
field of view. This contrasts with the 0.04 inches expected from the 
manufacturer's specification for the number of independent values in a line. 
Another estimate of the resolution may be obtained from a comparison of 
expected and observed noise levels. The expected noise level is calculated to 
be 2.15 Grinnell units (the scale used for temperature in Figure 3.3), while 
an analysis of the noise observed in one of the frames used in producing 
Figure 3.3 gives a standard deviation of 3.5 units. The expected noise value 
is based on the value of 100 independent samples in the horizontal field of 
view. If this factor is increased to a value consistent with 32 independent 
values per line (V(100/32)), the expected noise becomes 3.8 units, quite close 
to the observed 3.5 units. Thus, it appears that the horizontal resolution of 
the imager in this data was less than the manufacturer's suggested value by a 




factor of 3. The difference may be due to the employment of the electronic 
zoom feature of the imager to fit the video image into the full video frame, 
thereby enlarging the portion of the image occupied by the electronic 
footprint. Other possibilities are poor focusing in the data or vibration of 
the imager relative to the sample. Since the noise analysis and the step 
function analysis produce similar values' for system resolution, it seems 
appropriate to investigate the system resolution directly in a continuation of 
this work. 

In another comparison, the value of the vertical axis corresponding to 
the far-field stress used in the model is 18 of the Grinnell units. An 
expected value can be' calculated from handbook values of thermal expansion 
coefficient, heat capacity and density, using equations 3.1 and 3.2. The 
expected value is 18.8 of the Grinnell units. Given the uncertainty in 
determination of the model parameters, estimated at 2 Grinnell units, the 
agreement here is excellent . 

Finally, the figure of the experimental curve is more "rounded" on the 
top of the experimental data than on the bottom. One effect of this rounding 
is that the peak value of the observed data, corresponding to the point of 
maximum stress, misses the expected value by the greatest amount of any part 
of the image. This is of concern because it is this value, the greatest 
stress in the pattern, which is of greatest interest in much engineering 
stress analysis. In addition, the data lie consistently above the model in 
the region of the top curve outside of the rounded top. The roundedness of 
the top with associated the vertical offset further out are of some interest, 
as one of the possible explanations is that the sample was taken slightly 
beyond the elastic limit with the 5000 pounds of force applied. With the cross 
section of .5 square inches, the far field stress is 10,000 pounds per square 
inch (10 Kips) . The "stress factor" for the infinite geometry is 3.00, where 
stress factor is the standard engineering term for the ratio of maximum stress 
to the far-field stress. When the first standard correction for width is 
taken into account, the stress factor for our plate is 3.13, resulting in a 
stress of 31.3 Kips for our configuration. When the finite width is taken 
further into account, the stress factor in the center plane of the plate is 
slightly larger than this value. The handbook value for 0.2% yield for the 
sample material is 35 Kips. Given the closeness of the experimental stress to 
the 0.2% yield value, the character of the yield curve for Type 327 Stainless 
Steel, which is gradual, and the possibility that the experimental protocol 
may cause a little overshoot in the MTS machine, it is not unreasonable to 
consider that the high stress portions of the samplemay have yielded slightly. 
Measurements of the thickness of the inside of the hole revealed less than 
.001 inches of thinning, the limit of the measurement, but the expected yield 
would draw the material less than this amount. Because the sample had been 
subjected to more than 150 cycles of stresses at the 5000 pound level (but 
none at higherloads), the signal from initial yielding is not included in the 
analyzed data . The evidence for yielding would be present in the sample only 
as a residual stress. 

Under the experimental protocol, the sample was subjected to multiple 
cycles of stress and relaxation, the stress level being uniform for each 
cycle. The yielding would occur, and produce heating, during the first of the 
multiple cycles, but not on following cycles. Following the first cycle, the 
yielding would leave a pattern of residual stress, which would be sufficient, 
along with the work hardening, to halt the strain on the subsequent cycles 
just prior to the yield point. In conversations, with analysts, it appears 
that the strain pattern subsequent to the yielding would be different from 
that without yielding, and that the resulting alterations in the strain 
pattern would qualitatively resemble those in the observations. It is thus 
plausible that the rounding of the upper peak in the data is the result of 
residual stresses in the sample due to yielding during the experimental setup. 

The detection of residual stress or evidence of prior excedence of the 



elastic limit is of interest in many potential applications of non-destructive 
testing. It is recommended that the experiment reported here be continued by 
subjecting a sample similar to the one examined here to multiple cycles at 
increasing stress levels to see if the rounding pattern occurs subsequent to, 
but not prior to yielding and to examine its evolution as the load is 
increased to produce more severe yielding. 

3.3 Application to Composite Materials 

Samples of impacted composite material have been placed in the apparatus 
and examined both with' our imaging system and with an Ometron Spate 8000. The 
direct result is that the thermoelastic signal is much smaller than in steel 
samples and that direct contrast around the impact damage is not much larger 
than that in the region not subject to impact. On the other hand, when an 
examination of the region centered on the impact point is made, the 
thermoelastic signal does seem to change qualitatively in its character 
compared to an unimpacted region. Thus, the utility of thermoelasticity to 
detect impact damage in composite materials is neither immediately established 
nor clearly disproven. To examine the question further, it is necessary to 
look further into the nature of thermoelasticity in graphite-epoxy composite 
materials. 

At the most primitive level, considering only the bulk properties of the 
materials, many of the fabrication patterns for graphite-epoxy are designed to 
have a very small, essentially zero, coefficient of thermal expansion. This 
coefficient is in the numerator of the expression for the thermoelastic 
coefficient, so a zero value implies that there will be no thermoelastic 
temperature response to stress. This provides some insight into the small 
levels of thermoelastic signal observed in graphite-epoxy samples. The 
signals are not zero, however, and patterns of contrast are seen in 
experimental data. 

At the next level of complexity, a graphite-epoxy material may be 
considered to be a composite of two materials, graphite fibers embedded in an 
epoxy resin matrix. The graphite fibers are highly anisotropic, having very 
high strength and modulus of elasticity in their long direction. Their 
coefficient of thermal expansion in the along-fiber direction is actually 
slightly negative at room temperatures . In the cros3-f iber direction, their 
coefficient of thermal expansion is actually fairly large and positive. On 
heating, then, there is a net volume expansion of a graphite fiber in such a 
way that it becomes slightly shorter and a lot fatter. In contrast, the epoxy 
resin matrix is relatively isotropic, weak, and even compressible, with a 
fairly large coefficient of thermal expansion. If one considers the 
transmission of stress along the fiber direction in a composite material made 
of these two components, the majority of the load in a section is carried by 
normally stressed fibers, while the matrix transfers load between adjacent 
fibers through shear stresses at relatively low levels. The thermoelastic 
effect is associated in isotropic materials only with normal stresses. The 
thermoelastic effect from the fibers is thus zero, as their coefficient of 
thermal expansion in the load direction is zero, and that from the matrix is 
also zero, as it is loaded in shear. One can envision loading a composite 
material at an angle to its primary axis and producing normal stresses in the 
matrix material, which then give rise to thermoelastic signals. With this 
perspective, it may be that the generation of thermoelastic signals which are 
observed in composite samples is related to non-uniformities in the 
distribution of stress within the material which lead to patterns of normal 
stress in the matrix. Observation of such non-uniformities may be of interest 
in NDE and process control. 

Two further factors need to be considered before the observed patterns 
can be related to material flaws. The observed signals come primarily from 
temperature changes generated at the surface of the material, with the 



associated contrast transmitted to the surface layer by means of the stress 
distribution- The signals arise only secondarily from temperature anomalies 
generated within the bulk of the material and transmitted to the surface by 
thermal conduction. These two paths of transmission may be discriminated 
through the delay of the thermal conduction signal relative to the 
stress-transmitted signal. To the extent that the observed signal is produced 
to the dilitation of the surface layer, techniques other than 
thermoelasticity, involving special surface coatings which are brittle or 
change optical properties in response to strains, may provide superior NDE 
probes. The other factor' is that the theory presented here is based on 
isotropic materials while the graphite fiber component of the composite is 
notably anisotropic. In order to obtain an independent estimate of the 
thermoelastic temperature from the graphite component of the composite, the 
thermoelastic theory for non-isotropic materials must be found or derived and 
applied to graphite fibers. 

It is recommended that further research into the application of 
thermoelasticity to the NDE of graphite-epoxy composite materials take the 
path of describing the source of the observed signals. From the present, 
incomplete, perspective, it appears that the signals obtained from 
thermoelastic measurements arise from non-uniformities in the distribution of 
strain within a material rather than as characteristics of the normal state of 
the material subjected to cyclic stress. 
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Appendix 1 

Summary of Computer Programs 


Thermal Models 

POINTEMP.FOR - This program uses a Green's Function analysis to compute 
the surface temperature as a function to radius and time due to a 
Gaussian flux of heat absorbed on the surface of a thin plate over a 
finite length of time at a uniform rate. Zero heat loss is assumed at 
the front and back face of the plate. Physical and control variables are 
stored in two separate files COMTROL . FOR and PHYSICS. FOR. The function 
routine ERFC is called to produce values for complementary error 
functions used in the calculation. The resulting tables, stored in an 
array called POINTPLOT .DTA can be expressed in a form compatible with the 
plotting package IAP with program POINTPLOT. FOR;' Further calculations of 
estimates for spatial and temporal derivatives used for estimating 
diffusivity by the primitive equation method (see text) are done on the 
POINTOUT . DTA files with the program CALCK.FOR and placed in another file 
called DIFFEST. DTA. Physical parameters, temporal increments and limits, 
-pulse durations and radial scale, and radial increments and limits are 
all chosen by the user. 

LTEMP . FOR - This file contains one of the masters for program LINETEMP. 
This program uses Green's function analysis to produce the estimated 
thermal pattern at a given time due to a Gaussian laser source which is 
scanned in a straight line over the face of a semi-infinite (very thick) 
material. The version stored in LTEMP. FOR produces the estimate for an 
image with 480 lines each with 512 pixels. The version stored in 
SMLLTEMP.FOR uses an image of 120 lines and 128 pixels. The smaller 
version was found to run much faster than the larger one (see text) and 
is recommended for use. File references and parameter values are stored 
in files FILE. FOR and PARAM.FOR, with values obtained in subroutines 
GETFILES . FOR and GETP ARAMS . FOR . (For SMLLTEMP, GETPARAMS . FOR is replaced 
with LPARAM . FOR . ) Calculation boundaries are calculated in subroutine 
BOUND. FOR, and calculation of laser (heat pulse) intensity is done in 
subroutine PS. FOR. Pulse parameters that are user-specified include 
power, line position, line length, line scale width, line orientation, 
and pulse length. Material properties include thermal conductivity, heat 
capacity, and surface emissivity. Image parameters include pixel size 
and image time. The output is in a form that can be translated into the 
form used in the image processing package SADIE. 

TLINE .FOR is a subroutine that calculates the temperature in a thin plate 
of homogeneous material subsequent to an initial temperature which 
represents the cross-section of a linear temperature pattern extending 
through the thickness of the material. It uses an evaluation of Laplaces 
solution to the one-dimensional heat equation (see text) . The subroutine 
calls the function ERFC. FOR to evaluate the value for the complementary 
error function. 

LINEVFY.FOR is a program which calls TLINE and uses it in an 
initialization-prediction mode to compare with data from actual 
experiments. The routine contains two user-chosen parameters, the 
thermal diffusivity of the material and a heat loss factor, corresponding 



to the rate of fractional heat loss from the sample through processes of 
radiation and convection. Experimental data is retrieved from files 
written in the "9 line format" by subroutine GET_TEMP_DATA, supplied by 
W. Winfree. The program uses the first line of the 9 experimental lines 
as initialization and producesdata, model and difference lines in a form 
directly compatible with the plotting package IAP. In doing its 

-comparison, it removes a baseline value. for- each. data line baaed on the 

average of 50 values (out of 256) from each end, calculates the total 
heat at each time, calculates the total variance as a function of time, 
and passes the time data through to an IAP file. 

RATIO. FOR-- This program evaluated the thermal pattern due to a point 
source (see text) on a thin but finite plate and showed that the 
primitive equation diffusivity calculation produced no region of 
asymptotically correct values in the vicinity of the source. 

Thermoelasticity Programs 

MTS. FOR - This program, written by A. Dail, controls the experimental 
timing between the VAX computer, the Grinnell Image Processor and the 
LPA-11 laboratory interface to execute a thermoelastic experiment. Many 
of the experimental parameters (see text) are chosen by the user for each 
experiment. - It saves the resulting image file in a user-defined file, 
which can be changed into ~a "format compatible- with the image analysis 
-■ ------- package -SADIE with, program S ADD AT, written by J. Moore. 

MOD. FOR - This program calculates model values for thermoelastic 
temperature fields along the central transverse line across a plate with 
a hole in it (see text) . It filters the data with a boxcar filter to 
simulate the effect of a finite "footprint" for the active detector 
area. 

FIT. FOR - This program compares the calculations from a model to the data 
at the same set of position coordinates and calculates the variance of 
the residual . 

Miscellaneous and Utility Programs 

DISPLAY. FOR - This program, written by A. Dail, displays data in an image 
file format on the Grinnell Digitizer viedo output. 

TEMPDIS .FOR - This program, written by A. Dail, displays a user-chosen 
"window" of data in an image file format as a pair of hexidecimal digits 
on the screen of a terminal. It is used for low-level tracking of 
numbers through the video imaging process . 

ERFC . FOR - This function subroutine provides a fast calculation of 
complementary error function used in many thermal models. It is extended 
analytically to provide symmetric values for negative values of the 
argument as erfc(-x) = 2.0 - erfc(x). 

BIFR.FOR - This is a routine for the I Bessel Function (Modified Bessel 
Function of the First Kind) which was imported from the NASA Central 
Computer Facility and changed to operate on the VAX. 



Appendix 2 

Translation of "Infrared Thermography and Structural Mechanics" 

During the literature review, the paper "T616thermographie infrarouge et 
m^canique des structures" by B. Nayroles,_e£. al . was discovered describing 
French work which is closely related to our work both in approach and in 
purpose . —To- provide this information to the team within MC1S, the paper was 
translated into English as -part of-the effort under this contract. The 
translation follows . 




INFRARED THERMOGRAPHY AND STRUCTURAL MECHANICS 

by 

B. Nayroles, R. Bouc, H. Caumon and J. C. Chezeaux 
C.N.R.S Laboratory of Mechanics and Acoustics, Structural Radiation Group, 
31 Joseph-Aiguier Rd., BP 71, 13277, Marseille Cedex 9, France. 

and 

E. Giacometti 

National Society for Aerospace Industry, (S.N.I.A.S) Helicopter Service, 

Marignane, France 

Abstract 

Possible applications of recent IR movie systems are looked at in the 
field of structural mechanics . Two main phenomena are involved, on the one 
hand, the isentropic temperature variations which occur during a periodic 
evolution, on the other hand, the dissipation of mechanical energy. After 
some elementary recalls of thermodynamics, some details are given concerning 
the accuracies which are available when a suitable signal processing is 
undertaken or not. Finally the largest developments seem to be possible in 
the next years . 


1. INTRODUCTION 


For several years, infrared cameras, which provide visual images of 
temperature variations less than 0.2°C around ambient, have been commercially 
available. Whatever the method of detection, a photo-sensitive cell or a 
pyroelectric detector, these instruments deliver an analog picture signal and 
synchronization signals which allow reconstruction of an image obtained by 
scanning: the standard conf iguations provide information which can produce 
immediate visual images, if attached to displays, or optionally may be 
digitized. In fact, the most recent uses are instantaneous visual inspection 
of heat flow in buildings, overheating in electrical circuits, or other 
thermal anomalies, such as those due to subcutaneous tumors. It is remarkable 
the substantial industrial and medical market already open for these devices 
while their laboratory use remains quite limited. 

It is natural to inquire what applications these devices might have 
in the Mechanics of Solids, either in their present state or with improvements 
attainable through numerical processing. This is the type of question we 
propose to answer in this article along with presentation of results already 
obtained in our laboratory and those of S.N.I.A.S, with which we have been 
associated for a large portion of these studies. 

While the literature is generally sparce, we cite works of particular 
interest from fracture mechanics [1,2], non-destructive testing [3-5], 
vibration mechanics [6] and plastic deformation [7] . Finally, in an article 
to appear [8,9], an experimental study of the thermomechanical behavior of 
materials is considered using stroboscopy of IR thermography. First, we 
consider the source of heat caused by dissipation of mechanical energy: that 
due to internal friction present within the structure or that due to plastic 
deformation. Nevertheless, the reader should realize that the experiments 
described are only feasability demonstrations, a large number of physical 



parameters having not been measured or determined with the requisite 
preparation. In contrast/ a substantial experimental program is planned for 
the near future to compare experiment with theory for dissipation mechanisms 
at a crack tip (Nguyen Quoc Son [10]} 

We will see that the application of infrared thermography to structural 
mechanics is restricted because of the limited precision available in the 
standard commercially available models. On the other hand, digitization of 
the signal, summation of the images and the resulting improvement in 
signal-to-noise ratio permits considerations such that, with such accuracy as 
well as speed of operation, or ease of use, IR thermography should gain an 
important place in structural mechanics laboratories. 

2. PHENOMENA OF THERMOMECHANICAL STUDIES 
2.1 Review of the three-dimensional equations of thermomechanics 

We consider a three-dimensional solid which occupies a volume ft in space 
with its exterior surface denoted by dft. Under the hypothesis of 
inf ititesimal geometric displacements, one obtains the relation 

e ij = (1/2) (u if ) + u jti ) (1) 

where £ designates the deformation tensor and u the vector displacement. 

We make the general hypothesis that the various energy variables are 
distributed, i.e. they can be represented by volume densities as follows: p, 

mass density; (1/2) pu 2 , kinetic energy density, U, internal energy density; 

W, Helmholz Free Energy density; S, entropy density; D, density of power 
dissipation. 

Further, we suppose that at all points M within ft, the Helmholz Free 
Energy, W, is a function of the deformation, E, the absolute temperature, T, 
and "internal variables." The formalism used here is applicable, without 
particular limitations, provided the internal variables are the components of 
a variable, i], which can be expressed in a vector space of finite dimension. 
To illustrate, one may consider T} to be a tensor of the same type as E. This 
implies, in summary. 


W = W(x,E,T},T) . 

Furthermore, by the definition of W 

° = le* s = u = w + TS (2) 

where c denotes the stress tensor. 

Designate respectively by P[<ec anc * p cal the mechanical and thermal power 

supplied to the solid by its surroundings. 

The First Law of Thermodynamics is written 

J Q <« + puJidn- p mec + p CAL «> 

while the Second Law is written 

P CAL - J Q (TS - D)dflt, D£0 (4) 

Moreover, the mechanical energy balance is 



Substituting in (3) P CAL and P MEC by their expressions in (4) and (5) and 
using relations (2) to expand U, one obtains conventionally. 


D 



( 6 ) 


Next, denote by q the vector flux of heat within the solid, and assume 
the linear law of Fourier 


q = - k • grad T 


(7) 


where k is the tensor of heat conduction. Denote by <j> the volume density of 
any possible sources of heat within £2; we can express 


CAL 


n -q dQQ) + [ <{> 

a J ci 


where n designates the unit outward normal vector. The divergence theorem 
gives 

P CAL = ♦ - div <*> dfl (8) 

On replacing P CAL by its expression in (4), and seeing that the relation 
should apply to all sub-volumes of £2, one obtains the local relation 


<J> + div(k- grad T) = TS - D. 


Now, calculating S from the second of equations (2) produces 

( 2 2 2 

.g + sLJl .fj + t \ 
3t3e 3t3*1 m 2 ) 


On the right-hand side, the tensor quantities 

i . ■ 1 a!* m = _ alii 

3t3e ' 3T3n 

are the thermoelastic coefficients, while the scalar 


e.Tl 


= _ t 21h 

3t 2 


is the volumetric specific heat with E and TJ constant. 
With these notations, the heat equation is written 


(9) 


c T = div(k-grad T) - (1*6 + m>f|)T + + D. (10) 

E,T| 

Examination of the right-hand side reveals that temperature changes are 
due to three types of causes: 

1. Conduction, represented by the first term, which tends to make the 
field uniform 

2. Reversible effects, represented by (l-£ + m*^); when conduction, 

dissipation and internal sources can be neglected, the temperature varies 
with a linear combination of apparent and hidden deformations; the 
coefficients 1 and m are small, but the sensitivity of existing cameras 
permits observations of these variations. 


3. Heat sources, with which the dissipation of mechanical energy may be 
grouped. 

It is the reversible terms and the dissipation which we wish to evaluate 
from measurements of the temperature field made by an IR camera. 

With equation (10) is associated the equation of motion, 

pii + div o = f, (11) 

the geometric relation (1) and the equation of state, symbolically 


a = F (e,T) (12) 

where the functional F is in fact defined by a differential system depending 
on o, e, T\ and T. We note that the derivative of the first of relations (2) 
gives 

it - ,i3 > 

To illustrate the developments which follow, we present two examples 
which are simple as well as highly useful in practice 


Eaample 1: Linear Viscoelasticity 

We consider a substance whose mechanical behavior is represented by the 
following rheological model (Fig. 1) . The upper portion represents the 
thermal dilitation e°(T). Define the tensors of stiffness and the viscosity 
K 1 , K 2 and V, considered independent of temperature, and the specific heat at 
constant stress, C, which are relatively easily determined experimentally. 
The equations of mechanics are 

a = K 2 - [e - t) - e ° (T) ] = K 1 * T\ + V-i). 

From the first and the obvious expression 

D = i) -V -f| 

one infers that the form of the Helmholz Free Energy is 

w ( e - q - e°(T)] -K 2 ■[ e - q - e° (T) ] + + f <T) ( 14 ) 


where the function f remains to be determined. At zero strain, e is 
reduced to e°, and (10) reduces, assuming f = cT and grad T <= 0, to 



= -T1 


• T +CT 
dT 


since from (13) one infers 

1 = -m 


K 


del 

dT 


(15) 


it leads to 



. o ,o 

T d£_ . 2 d£_ 

dT K 


dT 


-Tf" 


and finally 

w = j [ e- q - e°(T) ] -k 2 -I e - -q - e° (T) ] + ^r\-K 1 -r\ + cT(l - log T) . ( i 6 ) 


Example 2. Standard Elastoplasticity 

The case of rheology is sufficiently similar: it suffices to replace K 1 

by 0 and V by a rigid element with plastic deformation t), more traditionally 
denoted £ p . One thus obtains 



D = a-£ p 

W = 1 [ e - e p - e°(T) ] -K 2 •[ e - e p - e°(T) J + cT(l - log T) . 

We return, now, to eq (10); the reversible terms evaluate in this case 
as before to 

— - T41^E--+- m-fj ) = J '4 t“‘K 2 ‘-( 

dT dT 

so that equation (10) is now written 

c ^ - div(k*grad T) = - T ■ b + D + (17) 

where in the first term, c is the volumetric specific heat with zero stress. 

This result is not limited to the preceding examples: it is easily shown 

that it extends to those solids for which the Helmholz Free Energy W is of the 
form 


W(£,r),T) = + f (T) , 


(18) 


that is, generally, to those which may be represented by a model of the 
preceding type with coefficients independent of temperature. One then has 


c 


c + fi£_° . 'filHi . de_? T 

e * T l dT de.de. dT 


2.2 Flat Plate Approximation 


An IR camera cannot, of course, measure the interior temperature of an 
object, but only its surface temperature. Therefore, in general, the 
information of significance to mechanics does not appear except by inference 
which follows the solution of an inverse problem. Nevertheless, the 
temperature may often be considered practically uniform through the thickness 
of a flat plate; one may then treat the situation as a two-dimensional solid, 
and an IR camera provides access to the temperature field in the plan domain. 

It is not our purpose here to discuss definitively the quality of the 
flat plate approximation in general, but simply to show that the equation for 
the thermal field in a homogeneous and isotropic flat plate subject to plane 
stess is of the form 



-a tr ( b ) + c 


(19) 


and to evaluate the order of magnitude of the different terms. 

We consider, in effect, such a plate (Fig. 3) of thickness e with its 
front face situated on the plane z=0; the tensor of conductivity and that of 
dilitation (de°/dT) are then scalars, and we can write 

a = Zte-r d = ^ (thermal diffusivity) 

We further consider our interest to be in small variations of temperature 
(a few degrees) around a reference temperature T 0 obtained, for example, when 
the plate is in thermomechanical equilibrium: stresses, strains and 

temperature are independent of time, and (17) is written as 



Since the coefficients a, c, and k vary little with temperature, and 
since the thickness is small, it is reasonable to assign to these 
coef f iciernts, at the point (x,y, z), the values which they have at temperature 



T o( x 'y/0)» thus, they become functions of x and y only. 
The difference, T' = T-T 0 is the solution of 




K dx 2 dz2' 


is ) + 


Now define 


e = iJ 


0 

T' <x, y, 2, t) dz 

-e 


( 20 ) 


and take the mean value of both sides of (20) 

The left-hand side produces the expression 

' 2/ ) r) 2 (i \ 1 I rim * 

(x,y,-e) 


3&- <*(£4 + a!£i) + j .j k &L± 

dt v 8x 2 By 2 ' ec L k 9z 



where the term within the square brackets is the flux of heat lost by the 
plate through the front and back faces. If T 1 is small, the square brackets 
are, in general, equal to a linear function of T’(x,y,-e) and T'(x,y,0) 
through linearization about T 0 of the laws of conduction, convection and 

radiation. Moreover, T' (x,y,-e) and T' (x,y,0) may be approximated by 0 such 
that the preceding expression is identical to the left-hand side of (19) . The 
constant T th has the dimension of time; in the case where the effects of 

convection and conduction may be neglected, this constant is given by the 
Stefan-Boltzmann law 

q = emiss x <T E T 4 

where q designates the radient flux, emiss the emissivity (equal to 1 for a 
black body) , and <r E is the constant 

O e « 5.67 x lO" 8 W/m 2 / 6 K 4 . 

On linearizing about T 0 , one obtains 



emiss. 


ec 


8 OgT 


( 21 ) 


which is, for a black body at 300K 

\h = ec x 12.25 sec. (22) 

For common materials, T th •= is on the order of hundreds of seconds per 

millimetre of thickness. When the plate is in contact with air, convection 
can reduce T th to around ec/10 or even to ec/50. Therefore, one may wish to 
perform experiments in a partial vacuum to eliminate convection. 

Under the plane stress hypothesis, the stress tensor o and that of strain 
£ are independent of z, while the components a 13 , a 23 and <J 33 are zero. The 

dissipation D is thus also independent of z, and finally, one obtains equation 
(19) . 


2.3 Some orders of magnitude. 


To fix these concepts, we assign numerical values to the various 
constants which appear in (19), for a black plate 1 mm in thickness and 
constituted of one or the other of the two quite different materials, titanium 
or plexiglas. The temperature T 0 is room temperature: 300K. To estimate the 
terms a tr (s) and D which occur in (19), one considers the two materials to 
be pulled within their linear domains, with the variations of stress C, to be 
simple tension or compression o m and -0 m , 0 m designating the elastic limit in 
tension. This permits the comparison of the two selected materials despite 
their great dissimilarity. We introduce tan $, $ designating the conventional 
loss angle, that, is the phase between strain and stress under harmonic motion 


of frequency f t . D is calculated under the hypothesis of such a motion at a 
single frequency f t remote from a "peak, of internal friction." One has 

D(t) = E -1 O m 2 2nf t tan $ sin 2 <2»cf x t) (23) 

and it is in fact the quantity 

-l 2 

80 = — n tan <j> (E: Young's Modulus) (24) 

c 

which is observed: this is the elevation in temperature caused by dissipation 
during one period in a volume element having zero heat loss. Likewise, the 
quantity 

50 = 2 aa m (25) 

represents the total temperature variation produced during one cycle by the 

"isentropic" term atr(ff). 

A 

One notes that the "isentropic" variations of temperature, that is d0 
are, during the course of a single cycle of magnitude 1000 times (titanium) or 

20 times (plexiglas) greater than the elevation of temperature 60 due to 
dissipation when “ the maximum stress is the elastic limit. 


Quantity 

Units 

Titanium 

Plexiglas 

c 

( j/m 3 ) /°K 

2.35xl0 6 

1.76xl0 6 

k 

(W/m) /°K 

21.75 

0.188 

d 

m 2 /s 

9.26xl0“ 6 

0 . 107xl0 _s 

3e°/3T 

°K _1 

8 . 9xl0~ s 

80xl0 -s 

a 

°K/pascal 

1 . 14xl0~ 9 

13 . 6xl0 -9 


hectobars 

35 

5 

E 

hectobars 

10550 

300 

tan <t> 

1 

0 . 5xl0 -3 

40x10*3 

80 

°K 

0. 78xl0- 3 

60x10*3 

80 

°K 

0.80 

1.35 

W e 

sec/mm 

192 

144 


These values apply at maximum stress; for smaller values, 80 varies with the 
maximuim stress and 80 with its square. 

3. OBSERVATIONS OF OSCILLATIONS OF A PLATE 

The preceding numerical values indicate sufficiently clear discrimination 
between the isentropic variations of temperature, which are observed on a time 
scale short with respect to l/f t and the average heating, which is observed on 

a time scale long with respect to T th , to the extent that a stationary regime 

is obtained for which the heat loss associated with temperature increase 
removes the heat generated by dissipation. 

We thus consider a flat plate subject to oscillating forces in its plane 
of frequency f t . If this is clearly lower than the frequency of the gravest 
plate mode, the mechanical force regime is attained in a few periods; if we 
assume further that the small temperature increase does not change the 
conditions for dissipation, by causing, for example, the onset of plastic 
deformation, one may further consider the right-hand side of (19), which we 



It is conventional 


shall call q, to be a periodic function with period f t -1 . 
to decompose this function as 


-l 


- - „ f't 

q = q + q/ q = ftj q(t)dt 


(26) 


One further finds 


q = -a tr(CT ) + £ - (27) 

Although the temperature 0 is not a priori a periodic function of time, 
one may still put 

-1 


t+f t 

rje($)d$ 


0 = e + e. e (t) = f 


(28) 


which separates the "slow" part 0 from the "fast" part 0 . 

One can show that (19) separates into two equations as follows 


30 

at 


- dV H 2 0 + 


‘th 



30 _ hV 2 £ + e = 

3F d H0 W 


= -a tr ( a ) + 


(29) 

(30) 


where V H denotes the bidimensional Laplacian 

2 2 
v/ __ 3 — x 3 

H “dx 2 dy 2 

This separation is valid under the conditions of linearity which we have 
assumed. 


Evolution of the mean temperature 0 

This is governed by equation (29). If one assumes that_ x t h is 
non-zero, then (29) certainly has one solution of constant 0 W . We note 
that if T th is zero, this existance may still be assured by virtue of the 
limits allowed on removal of heat supplied by D. Conventional inequality 
analysis then allows us to write 

II0OO-0 (t) II ^ exp [-—] lle^-efo) I (31) 

and an inequality of the same type (obtained using Gronwall's lemma) in the 
case where T th is zero and the temperature is fixed on the boundaries. 

Applying a spatial Fourier'*’ Transform F x y one can write equation (29) in 
the form. 


+ [d(f 
3t 


2 

x 


+ f v ) 




k th 




F 

C *x,y 


(D) 


and the same technique of inequality analysis produces, under this 


t In the case of a finite plate, one should, in all rigor, use a finite Fourier 
transform, that is, the Fourier transform of the product of 0 and a "bounding window" whose 
opening is occupied by the plate. One reaches the same conslusions. 



transformation 


I® 


«e,x,y 


- 0 x,y (t) II ^ exp ( - t [ d ( f x + fy ) H©..x. 

th 


© ( 0 ) 

x,y 


(32) 


generalizing (31) and showing that the limiting regime is attained very 
rapidly for the high spatial frequencies f x and f y . 


Evolution of the term with zero mean value 0 

It is interesting to compare the orders of magnitude of the terms 
directly. A spatial and temporal Fourier transform F x t of (30) thus 

produces 

|j> 2,t£ t + 4*d<f 2 + fy W)) + ir iy t (i|. 


(33) 


In the absence of diffusion, loss through the faces and dissipation, only 


the first term remains on each side of the equation, and the temporal 
variation of temperature proportional to that of stress is written as 

X 

^c,y,t = - a f t tr (f ( c) ) . 

z x,y,t 

—We evaluate the order of magnitude of the other terms as 
perturbations to these phenomena. It was previously shown that on the 
right-hand side, D was of the order of D since the temperature 
variations during a cycle are much smaller than the variations of the 
type 50- On the left-hand side, 1/T th is on the order of 10 , 

practically negligible in comparison with 2nf t .It is more interesting to 
calculate the temporal frequency f t , at which one can justify neglecting 
diffusion. 

Suppose that we observe an image of diameter 50 mm and that we want to 
distinguish details on the order of a half millimetre, that is, a very precise 
image . Thus 

■/fjP TT^ = 2xl0 3 iru 

This accomplished, if one wishes that the diffusion terms be about one 
percent of the leading terms, then 

f i> 87tdl0 6 = i 23200 Hz for titanium 
t V ; 1 270 Hz for plexiglas 

which are frequencies easily obtainable in the laboratory, and substantially 
below the lowest resonant frequencies, which are on the order of 100 KHz for a 
plate of titanium and ablut 31 KHz for one of plexiglas. This indicates the 
possibility of using the isentropic effect as a basis for a new technique of 
extensometry on thin objects which are subject to quasi-steady oscillations. 

Although the preceding estimates seem reassuring concerning the orders of 
magnitude of the various terms, a point yet remains to be settled: May one 

hope to obtain sufficient precision such that the theoretical possibilities 
which one glimpses correspond to potentially accessible techniques? It is 
this which we proceed to examine in the following paragraphs. 


4. VARIOUS EXAMPLES OF APPLICATION OF THERMOGRAPHY TO MECHANICS 

We wish first to describe an application for which thermomechanical 
coupling is not a cause: thermograms sometimes permit detection of certain 

anomalies in the thermal field for which the cause is geometrico-mechanical, 
as, for example, the presence of a crack. If this is open, thus essentially 



non-conducting to heat, the propagation of heat will conform to the shape of 
the obstacle, and so strongly perturb the temperature field.lt is this which 
one observes in the thermogram of figure 4, which represents a thin plate: a 

crack has been artificially introduced, and one clearly sees its effect on the 
thermal field since the temperature indicates a discontinuity. Applications 
of this type are already in the process of industrial development, and 
techniques of signal analysis and pattern recognition are expected to prove 
fruitful. 

The study of isentropic temperature variations has been pursued 
from the point of view of stroboscopy by Blanc and Giacometti [8]; the 
problem of precision undoubtedly requires a stroboscope at two points 
over one period ina fashion to eliminate the derivative of 8 . In 
effect, if one considers that the signal which gives 0 is swamped by 
noise of zero mean and standard deviation of about 0.1'C. one realizes 
that it requires about 100 measurements to bring the error to the 
vicinity of 1/100 ' C. The derivative of 0 is not more negligible, as 
that is the desired order of precision, even for a material with small 
dissipation. One sees that a device for rapid digitization of the 
signal is necessary. 

We expand further the applications to those which particularly involve 
the effects of dissipation. Various thermograms have already been published 
by [6J, who show plan view modes of vibration of diverse objects which are 
excited by ultrasonic transducers loudspeakers, plates in flexural vibration, 
etc. But the quantitative information supplied in these first published 
examples has remained very meager, no rapid digitization apparently having 
been employed. In contrast, we give in the following paragraph an example 
where, thanks to a relatively primitive digitizer, one may identify the 
sources of heat corresponding to plastic deformation in the vicinity of a 
crack tip. 


5. LOCATION OF THE PLASTIC DISSIPATION IN THE VICINITY OF THE TIP 

OF A FATIGUE CRACK 

We mounted a sample of titanium 60 mm wide and 0.5 mm thick on a tensile 
fatigue machine, and subjected it to cyclic oscillations at a frequency of 10 
Hz. A saw cut initiated a crack, the progression of which we followed with 
the aid of an AGA680 camera to obtain the photographs in Fig. 5 and 6, which 
show the field for the entire width. One distinguishes very quickly the 

heating 0 at the head of the crack, and the same photographs, which 
correspond respectively to the instants of minimum and maximum tension, 
also show the variation of q over the course of a cycle. The 
progression of the crack is very slow, on the order of 1 mm per 1000 
cycles; on the other hand, the spatial resolution (98 points per line, 70 
lines per image) is relatively small. In sum, on^ may consider the 
observed regime to be stationary and assume that q satisfies 

-d 7 H 2 e + -f— = E (34) 

H T th c * 

A 

One can thus evaluate the temperature 0 and eliminate 8 by an 
appropriate choice of sampling frequency, and the correct value for 
0(Xj.,yj) is, in fact, the arithmatic mean of N measurements of 0(x i ,yj,t) 
obtained at the point (x A ,yj ) for which the phase values are uniformly 
distributed between 0 and 2 ji. Here, N equals 100. The field of 
temperature 0 thus obtained is shown in figures 7 and 8. 



( 35 ) 



( 36 ) 


(see (11] for details.) 

The size of the significant extent of <{> provides the order of 
magnitude of the spatial resolution obtainable by this estimate. More 
precisely, these result from a low-pass filter ofVfl0 obs , as we see a 
little further on. 

By a classical property of convolutions, one has in consequence 
- - -- -7je es (x,-y) (X',y*) x e^U-xSy-y’Jdx’dy' 

which is descretized to the form 


V^8 es (x,y) =^> V* <)>(i'Ax, j'Ay) X e ob ( s (i-i')Ax, (j-j')Ay), (37) 

-Iii'SI 

-JSj'SJ 


the formulation leading to the results seen in Figure 9. 

2 — 

The estimate of (37) differs from the actual V H 0 by an error which 
arises from three sources. 


(a) Error from folding the spectrum (in English, aliasing) . It is 
well-known that sampling an actual signal leads to a systematic error if 
the sampling frequency is not at least two times greater than the highest 
frequency contained in the signal. _This type of error is much reduced 
here, for the field of temperature 0 has previously been concentrated in 
the lower spatial frequencies. Furthermore, the signal has been filtered 
by the digitizer. 

(b) Errors in filtering. The Fourier transform of (35) gives 

F l 7 H0est } < f x' f y> = F (^H0 obs Hfx,f y ) (38) 


As a simple example*, suppose that we have selected <I> = F(<J>) as the product of 
Hamming windows in f x and f y . Let be the one-dimensional window 


*H<P> 


{ 0.56 + 0.44 cos 2lp for 3) sfc- 

2 , 2 

0 for |p | 


(39) 


Accordingly, one takes f x e and f ye to be the spatial sampling 
frequencies along x and y. 

® (f xf f y) = (** /fx,e) / *y, e) 


The inverse transform to <I> is easily calculated, along with its 
Laplacian. One establishes that it is sufficient to limit, in (39), I=J=6, 
which defines the effective volume of calculation for the convolution. 
Elsewhere, one sees that (38) and (39) show the filtering error. If, for 
example, one accepts no more than 10% error of this type, one should be 
interested in only those spatial frequencies f x and f y such that 


(0.56 + 0.44 cos 2n 


fx 

fx.t 


) (0.56 + 0.44 cos 2X 




) 2; 0.9 


which gives, if one chooses the same limit of P max for the two axes f x /f x , e and 
*y/ f y,e : Pmax' 5 ®* 0 ®’ 

In past tests, the spatial frequency in x and in y has been on the order 
of 1.6 mm -1 , such that only details of diameter on the order of 
(1/2) (0.08 x 16) -1 «® 4 mm can be observed with a filtering error less than 10%. 
In optical terms, the resolution, for an error of 10%, is approximately a pair 
of lines every eight millimetres, which is weak. 


*In reality, one uses a window of revolution where the "meridian" is a Hamming window 



{ c ) Errors from noise. Consider each -temperature measurement to be 
associated with a random error with zero mean value and a standard deviation 
60 which applies to every measured point. Consider further that the errors 
between two different measured points are independent random variables. 


. 2 — • 

Then the estimate 7 U ft obtained with formula (37) is associated 

, n v est 

with a random error of zero mean and standard deviation 


- f i i/2 

5<V H 2 e) = Z[ 7 H<l>(i’Ax,j-Ay)] 2 | 

In our previous experience, as has been noted above, N=100 while 60 
is the raw resolution in temperature which, as was noted at the 
beginning of the study, is about 0.2 C. This gives for (5 a standard 
error of 0.02 C and, after application of (40), a standard error for 7^0 
of 


8<7 2 S> 


0 . 038 °C/mm 2 


(41) 


One may then write (34) below as 


- v hB + 


0 

<*th 


a 

k 


(42) 


“ • Even if natural convection considerably lowers x th with respect to 
values calculated for radiation from the faces, t th remains greater than 

0.2 sec for the cases studied. This limits the standard error 6(0 /dT th ) 
to 0.01°C/mm 2 . In all, the standard error for the left-hand side of 
(42) is less than 0.05°C/mm 2 . 


To estimate the quality of this approximation, suppose that, in the 
region under consideration of about 4 mm in diameter, the plastic deformation 



is, on the average, 10 3 per period. For one case studied at 10 Hz, that 
gives, for mean dissipated power 

D = 10 X 35 X 10 7 x 10" 3 W/m3 

and the right-hand side of (42) equals , , T 

^ = 0.161 ’c/mnf 


so that the standard error for the measurement of the right-hand side is on 
the order of 30%. 

Such a result is still very interesting in comparison with the 
possibilities of traditional techniques, as, for example, the diffraction %• 
But it needs to be considerably improved. 

1. By an optimal choice of optics, permitting the observer to concentrate 
on the zone of interest, thereby considerably reducing the filtering errors. 

2. With a rapid digitizer, to augment simultaneously the sampling 
frequency and the number N of obtainable measurements . 

3. To perform experiments under partial vacuum in order to 
eliminate convection losses and concomitantly the uncertainty in the 
terra (5 /dt th ), very difficult to control. 

One hopes, thus, instead of an analog experiment described here, to bring 
the standard error to a few percent for a spatial resolution on the order of a 
millimetre . 



Figure Captions 


Fig 1. Linear viscoelastic solid with dilitation. 

Fig 2. Ideal elasto-plastic solid with thermal dilitation e 0 . 

Fig 3. Thin plate "observed” on face B. 

Fig 4 Thermogram of a cracked plate 

Fig 5. Thermogram of the progression of a fatigue crack. IR stroboscopy at 
minimum stress. 

Fig 6. Thermogram of the progression of a fatigue crack. IR stroboscopy at 
maximum stress. 

Fig 7. Progression of a fatigue crack. Image is digitized and reassembled 
following noise removal . . 

Fig 8. (a) Thermal profile (noisy) along the crack axis. (b) The same profile 

.-averaged over ten frames, .(c) The same profile averaged over 50 frames, 
(d) The same profile averaged over 100 frames. 

Fig 9. Source of heat at the tip of a fatigue crack obtained by numerical 
filtering of the temperature data in Fig 7. 
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Fig. 5. Thcrmogramme dc progression dune fissure de fatigue. Stroboscopic IR sur I'effort minimal. 





iMSti 








i | | WMWI 




4 «; 


i' **$.■* *■«. ; 




’■ f 4 ~* , .|v | '.*• 

>$x% Wjj 


Fig. 7. Progression d'une fissure de fatigue. Image numerisee et reconstitute apres debruitage. 






